Data_Parame = readtable('parametros');
group=1;
while group<9
    if group==1
    Bootstrap=readtable('bootstrap');
    else
        file=['bootstrap',num2str(group),'.xlsx'];
        Bootstrap=readtable(file);
    end
    iter=1;
while iter<101
sigma=Data_Parame{group,1};
theta=Data_Parame{group,2};
rho=0.06;
w_w=1000*Data_Parame{group,3};
w_m=1000*Data_Parame{group,4};
r=Data_Parame{group,5};
kw_theta=(w_w/r)^(1/(1-theta));
mw_rho=(w_w/w_m)*((1+(w_w/r)^(theta/(1-theta)))^((theta-rho)/theta));
w_w=1.3*w_w;
w_m=1.3*w_m;
r=0.04;
RDKW=Bootstrap{iter,1};
RDKM=Bootstrap{iter,2};
frac_5=Data_Parame{group,6};
frac_10=Data_Parame{group,7};
frac_15=Data_Parame{group,8};
Mccrary=Data_Parame{group,9};

if RDKM<0 && RDKM/RDKW<0.3
    clear ratKM_jump

rho_grid=linspace(-6.01,-0.01 ,7);
theta_grid=linspace(0.01 ,0.31,7);

k=1;
while k<8;
    l=1;
    while l<8;
    rho=rho_grid(k);
    theta=theta_grid(l);
    lambda=(kw_theta)/(w_w/r+kw_theta);
    delta=mw_rho/(mw_rho+(w_w/w_m*(lambda^((1-rho)/(1-theta)))*((lambda^(1/(1-theta))+(((1-lambda)^(1/(1-theta)))*((w_w/r)^(theta/(1-theta)))))^((theta-rho)/(theta*(1-rho))))));
    ratwm=(w_w/w_m)^(rho/(1-rho));
    ratwr=(w_w/r)^(theta/(1-theta));
    ratwr2=(w_w/r)^(1/(1-theta));
    ratrm=(r/w_m)^(rho/(1-rho));
    ratrm2=(r/w_m)^(1/(1-rho));
    ratwm2=(w_w/w_m)^(1/(1-rho));
    l1=lambda^(1/(1-theta));
    l2=(1-lambda)^(1/(1-theta));
    l3=delta^(1/(1-rho));
    l4=(1-delta)^(1/(1-rho));

    alpha_5=(5^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_10=(10^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_15=(15^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

    flag=0;
    flag2=0;
    while alpha_5<1
        alpha_5=alpha_5*10;
        alpha_10=alpha_10*10;
        alpha_15=alpha_15*10;
        flag=flag+1;
    end
    while alpha_15>100
        alpha_5=alpha_5/10;
        alpha_10=alpha_10/10;
        alpha_15=alpha_15/10;
        flag2=flag2+1;
    end
     
    adjust=(10^flag)*(10^(-flag2));
    
    alpha_lbar=adjust*(19^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_20=adjust*(20^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_30=adjust*(30^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

       alpha_range=linspace(alpha_lbar, alpha_30, 1000);



    e1=(sigma-theta)/(theta*(1-sigma));
    e2=(rho*(theta-1))/(theta*(1-rho));
    e3=(sigma-rho)/(rho*(1-sigma));
    e4=(1-theta)/theta;
    e5=(1-rho)/rho;
    e6=1/(1-sigma);
    e7=(theta-rho)/(theta*(1-rho));

    W_lbar=(((1/adjust)*alpha_lbar*sigma*(l1^(1-sigma))/w_w)^e6)*((1-delta)^(sigma/(rho*(1-sigma))))*((l1+(l2*ratwr))^e1)*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^e2)))^e3);

    K_init=19*ratwr2*l2/l1;
    M_init=l3/l4*19*ratwm2*(1/l1)*((l1+(l2*ratwr))^e7);

    options = optimoptions('fsolve', ...
        'StepTolerance', 1e-15, ...  % Tolerance for the step size
        'FunctionTolerance', 1e-15, ...  % Tolerance for the function value
        'OptimalityTolerance', 1e-15, ...  % Tolerance for first-order optimality
        'MaxIterations', 10^10, ...  % Maximum number of iterations
        'MaxFunctionEvaluations', 10^10)  % Maximum number of function evaluations

    i=1;
    K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K_init, options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));

i=2;
        while i<1001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
            i=i+1;
        end

    %tau_grid=linspace(0.000001*w_w, 0.05*w_w,1000);
tau_grid=linspace(0,1500,1000);
    j=1;
    while j<1001
        tau=tau_grid(j);
        ratwtm=((w_w+tau)/w_m)^(rho/(1-rho));
        ratwtr=((w_w+tau)/r)^(theta/(1-theta));
        ratwtr2=((w_w+tau)/r)^(1/(1-theta));
        ratwtm2=((w_w+tau)/w_m)^(1/(1-rho));
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*((1/adjust*alpha_range*sigma*(l1^(1-sigma))/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            alpha_40=adjust*(40^(1-sigma))*(w_w/(sigma*l1^((1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
            alpha_range(1001:2000)=linspace(alpha_30, alpha_40, 1000);
            i=1001;
            while i<2001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
           i=i+1;
            end
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*l1*((1/adjust*alpha_range*sigma/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            I2=1001;
        end
        end
        if I2==2000
            I2=1999;
        end
        
        W_next=floor((1-delta)^(sigma/(rho*(1-sigma)))*l1*(((1/adjust)*alpha_range(I2+1)*sigma/w_w)^(1/(1-sigma)))*((l1+(l2*ratwr))^((sigma-theta)/(theta*(1-sigma))))*((1+(ratwm*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((sigma-rho)/(rho*(1-sigma)))));
alpha_next=adjust*(W_next^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
alpha_nextplus1=adjust*((W_next+1)^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

        ratKW_left(j)=mean(ratKW(I2:I2)');
        ratKW_right(j)=(((1-lambda)*(w_w+tau))/(lambda*r))^(1/(1-theta));
        ratKW_discon(j)=log(ratKW_right(j))-log(ratKW_left(j));
        %Aqui es donde tengo mas dudas
        ratKM_left(j)=mean(ratKM(I2:I2)');
        ratKM_right(j)=l4/l3*l2*(ratwtr2/ratwtm2)*((l1+(l2*ratwtr))^(-e7));
        ratKM_discon(j)=log(ratKM_right(j))-log(ratKM_left(j));
        j=j+1;
    end

    %diff_bunching=abs(bunching_density-Mccrary);
    diff_bunching=abs(ratKW_discon-RDKW);
    [D3,I3]=min(diff_bunching);
    tau_final(k,l)=tau_grid(I3);
    if D3<0.01
    ratKM_jump(k,l)=ratKM_discon(I3);
    ratKW_jump(k,l)=ratKW_discon(I3);
    bunch(k,l)=bunching_density(I3);
    else
    ratKM_jump(k,l)=9999;
    ratKW_jump(k,l)=9999;
    bunch(k,l)=9999;
    %ratKM_jump(k,l)=ratKM_discon(I3);
    %ratKW_jump(k,l)=ratKW_discon(I3);
    %bunch(k,l)=bunching_density(I3);
    end        
    clear pi_rest
    l=l+1;
    end
    k=k+1;
end

diff_KMjump=(ratKM_jump-RDKM).^2;
[D4,I4]=min(diff_KMjump);
[D5,I5]=min(D4');

rho_bygroup(iter)=rho_grid(I4(I5));
theta_bygroup(iter)=theta_grid(I5);
tau_bygroup(iter)=tau_final(I4(I5),I5);
bunching(iter)=bunch(I4(I5),I5);
KM_jump(iter)=ratKM_jump(I4(I5),I5);
KW_jump(iter)=ratKW_jump(I4(I5),I5);
end

if RDKM<0 && RDKM/RDKW>0.3 && RDKM/RDKW<0.6
    clear ratKM_jump

rho_grid=linspace(-15.01,-4.01 ,10);
theta_grid=linspace(0.01 ,0.31,7);

k=1;
while k<11;
    l=1;
    while l<8;
    rho=rho_grid(k);
    theta=theta_grid(l);
    lambda=(kw_theta)/(w_w/r+kw_theta);
    delta=mw_rho/(mw_rho+(w_w/w_m*(lambda^((1-rho)/(1-theta)))*((lambda^(1/(1-theta))+(((1-lambda)^(1/(1-theta)))*((w_w/r)^(theta/(1-theta)))))^((theta-rho)/(theta*(1-rho))))));
    ratwm=(w_w/w_m)^(rho/(1-rho));
    ratwr=(w_w/r)^(theta/(1-theta));
    ratwr2=(w_w/r)^(1/(1-theta));
    ratrm=(r/w_m)^(rho/(1-rho));
    ratrm2=(r/w_m)^(1/(1-rho));
    ratwm2=(w_w/w_m)^(1/(1-rho));
    l1=lambda^(1/(1-theta));
    l2=(1-lambda)^(1/(1-theta));
    l3=delta^(1/(1-rho));
    l4=(1-delta)^(1/(1-rho));

    alpha_5=(5^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_10=(10^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_15=(15^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

    flag=0;
    flag2=0;
    while alpha_5<1
        alpha_5=alpha_5*10;
        alpha_10=alpha_10*10;
        alpha_15=alpha_15*10;
        flag=flag+1;
    end
    while alpha_15>100
        alpha_5=alpha_5/10;
        alpha_10=alpha_10/10;
        alpha_15=alpha_15/10;
        flag2=flag2+1;
    end

     
    adjust=(10^flag)*(10^(-flag2));
    
    alpha_lbar=adjust*(19^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_20=adjust*(20^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_30=adjust*(30^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

       alpha_range=linspace(alpha_lbar, alpha_30, 1000);



    e1=(sigma-theta)/(theta*(1-sigma));
    e2=(rho*(theta-1))/(theta*(1-rho));
    e3=(sigma-rho)/(rho*(1-sigma));
    e4=(1-theta)/theta;
    e5=(1-rho)/rho;
    e6=1/(1-sigma);
    e7=(theta-rho)/(theta*(1-rho));

    W_lbar=(((1/adjust)*alpha_lbar*sigma*(l1^(1-sigma))/w_w)^e6)*((1-delta)^(sigma/(rho*(1-sigma))))*((l1+(l2*ratwr))^e1)*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^e2)))^e3);

    K_init=19*ratwr2*l2/l1;
    M_init=l3/l4*19*ratwm2*(1/l1)*((l1+(l2*ratwr))^e7);

    options = optimoptions('fsolve', ...
        'StepTolerance', 1e-15, ...  % Tolerance for the step size
        'FunctionTolerance', 1e-15, ...  % Tolerance for the function value
        'OptimalityTolerance', 1e-15, ...  % Tolerance for first-order optimality
        'MaxIterations', 10^10, ...  % Maximum number of iterations
        'MaxFunctionEvaluations', 10^10)  % Maximum number of function evaluations

    i=1;
    K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K_init, options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));

i=2;
        while i<1001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
            i=i+1;
        end

    %tau_grid=linspace(0.000001*w_w, 0.05*w_w,1000);
tau_grid=linspace(0,1500,1000);
    j=1;
    while j<1001
        tau=tau_grid(j);
        ratwtm=((w_w+tau)/w_m)^(rho/(1-rho));
        ratwtr=((w_w+tau)/r)^(theta/(1-theta));
        ratwtr2=((w_w+tau)/r)^(1/(1-theta));
        ratwtm2=((w_w+tau)/w_m)^(1/(1-rho));
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*((1/adjust*alpha_range*sigma*(l1^(1-sigma))/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            alpha_40=adjust*(40^(1-sigma))*(w_w/(sigma*l1^((1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
            alpha_range(1001:2000)=linspace(alpha_30, alpha_40, 1000);
            i=1001;
            while i<2001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
           i=i+1;
            end
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*l1*((1/adjust*alpha_range*sigma/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            I2=1001;
        end
        end
        if I2==2000
            I2=1999;
        end
        
        W_next=floor((1-delta)^(sigma/(rho*(1-sigma)))*l1*(((1/adjust)*alpha_range(I2+1)*sigma/w_w)^(1/(1-sigma)))*((l1+(l2*ratwr))^((sigma-theta)/(theta*(1-sigma))))*((1+(ratwm*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((sigma-rho)/(rho*(1-sigma)))));
alpha_next=adjust*(W_next^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
alpha_nextplus1=adjust*((W_next+1)^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

        ratKW_left(j)=mean(ratKW(I2:I2)');
        ratKW_right(j)=(((1-lambda)*(w_w+tau))/(lambda*r))^(1/(1-theta));
        ratKW_discon(j)=log(ratKW_right(j))-log(ratKW_left(j));
        %Aqui es donde tengo mas dudas
        ratKM_left(j)=mean(ratKM(I2:I2)');
        ratKM_right(j)=l4/l3*l2*(ratwtr2/ratwtm2)*((l1+(l2*ratwtr))^(-e7));
        ratKM_discon(j)=log(ratKM_right(j))-log(ratKM_left(j));
        j=j+1;
    end

    %diff_bunching=abs(bunching_density-Mccrary);
    diff_bunching=abs(ratKW_discon-RDKW);
    [D3,I3]=min(diff_bunching);
    tau_final(k,l)=tau_grid(I3);
    if D3<0.01
    ratKM_jump(k,l)=ratKM_discon(I3);
    ratKW_jump(k,l)=ratKW_discon(I3);
    bunch(k,l)=bunching_density(I3);
    else
    ratKM_jump(k,l)=9999;
    ratKW_jump(k,l)=9999;
    bunch(k,l)=9999;
    %ratKM_jump(k,l)=ratKM_discon(I3);
    %ratKW_jump(k,l)=ratKW_discon(I3);
    %bunch(k,l)=bunching_density(I3);
    end        
    clear pi_rest
    l=l+1;
    end
    k=k+1;
end

diff_KMjump=(ratKM_jump-RDKM).^2;
[D4,I4]=min(diff_KMjump);
[D5,I5]=min(D4');

rho_bygroup(iter)=rho_grid(I4(I5));
theta_bygroup(iter)=theta_grid(I5);
tau_bygroup(iter)=tau_final(I4(I5),I5);
bunching(iter)=bunch(I4(I5),I5);
KM_jump(iter)=ratKM_jump(I4(I5),I5);
KW_jump(iter)=ratKW_jump(I4(I5),I5);
end

if RDKM<0 && RDKM/RDKW>0.6
    clear ratKM_jump

rho_grid=linspace(-25.01,-10.01 ,16);
theta_grid=linspace(0.01 ,0.31,7);

k=1;
while k<17;
    l=1;
    while l<8;
    rho=rho_grid(k);
    theta=theta_grid(l);
    lambda=(kw_theta)/(w_w/r+kw_theta);
    delta=mw_rho/(mw_rho+(w_w/w_m*(lambda^((1-rho)/(1-theta)))*((lambda^(1/(1-theta))+(((1-lambda)^(1/(1-theta)))*((w_w/r)^(theta/(1-theta)))))^((theta-rho)/(theta*(1-rho))))));
    ratwm=(w_w/w_m)^(rho/(1-rho));
    ratwr=(w_w/r)^(theta/(1-theta));
    ratwr2=(w_w/r)^(1/(1-theta));
    ratrm=(r/w_m)^(rho/(1-rho));
    ratrm2=(r/w_m)^(1/(1-rho));
    ratwm2=(w_w/w_m)^(1/(1-rho));
    l1=lambda^(1/(1-theta));
    l2=(1-lambda)^(1/(1-theta));
    l3=delta^(1/(1-rho));
    l4=(1-delta)^(1/(1-rho));

    alpha_5=(5^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_10=(10^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_15=(15^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

    flag=0;
    flag2=0;
    while alpha_5<1
        alpha_5=alpha_5*10;
        alpha_10=alpha_10*10;
        alpha_15=alpha_15*10;
        flag=flag+1;
    end
    while alpha_15>100
        alpha_5=alpha_5/10;
        alpha_10=alpha_10/10;
        alpha_15=alpha_15/10;
        flag2=flag2+1;
    end

         
    adjust=(10^flag)*(10^(-flag2));
    
    alpha_lbar=adjust*(19^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_20=adjust*(20^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_30=adjust*(30^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

       alpha_range=linspace(alpha_lbar, alpha_30, 1000);



    e1=(sigma-theta)/(theta*(1-sigma));
    e2=(rho*(theta-1))/(theta*(1-rho));
    e3=(sigma-rho)/(rho*(1-sigma));
    e4=(1-theta)/theta;
    e5=(1-rho)/rho;
    e6=1/(1-sigma);
    e7=(theta-rho)/(theta*(1-rho));

    W_lbar=(((1/adjust)*alpha_lbar*sigma*(l1^(1-sigma))/w_w)^e6)*((1-delta)^(sigma/(rho*(1-sigma))))*((l1+(l2*ratwr))^e1)*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^e2)))^e3);

    K_init=19*ratwr2*l2/l1;
    M_init=l3/l4*19*ratwm2*(1/l1)*((l1+(l2*ratwr))^e7);

    options = optimoptions('fsolve', ...
        'StepTolerance', 1e-15, ...  % Tolerance for the step size
        'FunctionTolerance', 1e-15, ...  % Tolerance for the function value
        'OptimalityTolerance', 1e-15, ...  % Tolerance for first-order optimality
        'MaxIterations', 10^10, ...  % Maximum number of iterations
        'MaxFunctionEvaluations', 10^10)  % Maximum number of function evaluations

    i=1;
    K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K_init, options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));

i=2;
        while i<1001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
            i=i+1;
        end

    %tau_grid=linspace(0.000001*w_w, 0.05*w_w,1000);
tau_grid=linspace(0,1500,1000);
    j=1;
    while j<1001
        tau=tau_grid(j);
        ratwtm=((w_w+tau)/w_m)^(rho/(1-rho));
        ratwtr=((w_w+tau)/r)^(theta/(1-theta));
        ratwtr2=((w_w+tau)/r)^(1/(1-theta));
        ratwtm2=((w_w+tau)/w_m)^(1/(1-rho));
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*((1/adjust*alpha_range*sigma*(l1^(1-sigma))/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            alpha_40=adjust*(40^(1-sigma))*(w_w/(sigma*l1^((1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
            alpha_range(1001:2000)=linspace(alpha_30, alpha_40, 1000);
            i=1001;
            while i<2001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
           i=i+1;
            end
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*l1*((1/adjust*alpha_range*sigma/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            I2=1001;
        end
        end
        if I2==2000
            I2=1999;
        end
        
        W_next=floor((1-delta)^(sigma/(rho*(1-sigma)))*l1*(((1/adjust)*alpha_range(I2+1)*sigma/w_w)^(1/(1-sigma)))*((l1+(l2*ratwr))^((sigma-theta)/(theta*(1-sigma))))*((1+(ratwm*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((sigma-rho)/(rho*(1-sigma)))));
alpha_next=adjust*(W_next^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
alpha_nextplus1=adjust*((W_next+1)^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

        ratKW_left(j)=mean(ratKW(I2:I2)');
        ratKW_right(j)=(((1-lambda)*(w_w+tau))/(lambda*r))^(1/(1-theta));
        ratKW_discon(j)=log(ratKW_right(j))-log(ratKW_left(j));
        %Aqui es donde tengo mas dudas
        ratKM_left(j)=mean(ratKM(I2:I2)');
        ratKM_right(j)=l4/l3*l2*(ratwtr2/ratwtm2)*((l1+(l2*ratwtr))^(-e7));
        ratKM_discon(j)=log(ratKM_right(j))-log(ratKM_left(j));
        j=j+1;
    end

    %diff_bunching=abs(bunching_density-Mccrary);
    diff_bunching=abs(ratKW_discon-RDKW);
    [D3,I3]=min(diff_bunching);
    tau_final(k,l)=tau_grid(I3);
    if D3<0.01
    ratKM_jump(k,l)=ratKM_discon(I3);
    ratKW_jump(k,l)=ratKW_discon(I3);
    bunch(k,l)=bunching_density(I3);
    else
    ratKM_jump(k,l)=9999;
    ratKW_jump(k,l)=9999;
    bunch(k,l)=9999;
    %ratKM_jump(k,l)=ratKM_discon(I3);
    %ratKW_jump(k,l)=ratKW_discon(I3);
    %bunch(k,l)=bunching_density(I3);
    end        
    clear pi_rest
    l=l+1;
    end
    k=k+1;
end

diff_KMjump=(ratKM_jump-RDKM).^2;
[D4,I4]=min(diff_KMjump);
[D5,I5]=min(D4');

rho_bygroup(iter)=rho_grid(I4(I5));
theta_bygroup(iter)=theta_grid(I5);
tau_bygroup(iter)=tau_final(I4(I5),I5);
bunching(iter)=bunch(I4(I5),I5);
KM_jump(iter)=ratKM_jump(I4(I5),I5);
KW_jump(iter)=ratKW_jump(I4(I5),I5);
end

if RDKM>0
    clear ratKM_jump

rho_grid=linspace(0.01,0.11,11);
theta_grid=linspace(0.01,0.31,11);

k=1;
while k<12;
    l=1;
    while l<12;
    rho=rho_grid(k);
    theta=theta_grid(l);
    lambda=(kw_theta)/(w_w/r+kw_theta);
    delta=mw_rho/(mw_rho+(w_w/w_m*(lambda^((1-rho)/(1-theta)))*((lambda^(1/(1-theta))+(((1-lambda)^(1/(1-theta)))*((w_w/r)^(theta/(1-theta)))))^((theta-rho)/(theta*(1-rho))))));
    ratwm=(w_w/w_m)^(rho/(1-rho));
    ratwr=(w_w/r)^(theta/(1-theta));
    ratwr2=(w_w/r)^(1/(1-theta));
    ratrm=(r/w_m)^(rho/(1-rho));
    ratrm2=(r/w_m)^(1/(1-rho));
    ratwm2=(w_w/w_m)^(1/(1-rho));
    l1=lambda^(1/(1-theta));
    l2=(1-lambda)^(1/(1-theta));
    l3=delta^(1/(1-rho));
    l4=(1-delta)^(1/(1-rho));

    alpha_5=(5^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_10=(10^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_15=(15^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

    flag=0;
    flag2=0;
    while alpha_5<1
        alpha_5=alpha_5*10;
        alpha_10=alpha_10*10;
        alpha_15=alpha_15*10;
        flag=flag+1;
    end
    while alpha_15>100
        alpha_5=alpha_5/10;
        alpha_10=alpha_10/10;
        alpha_15=alpha_15/10;
        flag2=flag2+1;
    end

    
    adjust=(10^flag)*(10^(-flag2));
    
    alpha_lbar=adjust*(19^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_20=adjust*(20^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
    alpha_30=adjust*(30^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

       alpha_range=linspace(alpha_lbar, alpha_30, 1000);



    e1=(sigma-theta)/(theta*(1-sigma));
    e2=(rho*(theta-1))/(theta*(1-rho));
    e3=(sigma-rho)/(rho*(1-sigma));
    e4=(1-theta)/theta;
    e5=(1-rho)/rho;
    e6=1/(1-sigma);
    e7=(theta-rho)/(theta*(1-rho));

    W_lbar=(((1/adjust)*alpha_lbar*sigma*(l1^(1-sigma))/w_w)^e6)*((1-delta)^(sigma/(rho*(1-sigma))))*((l1+(l2*ratwr))^e1)*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^e2)))^e3);

    K_init=19*ratwr2*l2/l1;
    M_init=l3/l4*19*ratwm2*(1/l1)*((l1+(l2*ratwr))^e7);

    options = optimoptions('fsolve', ...
        'StepTolerance', 1e-15, ...  % Tolerance for the step size
        'FunctionTolerance', 1e-15, ...  % Tolerance for the function value
        'OptimalityTolerance', 1e-15, ...  % Tolerance for first-order optimality
        'MaxIterations', 10^10, ...  % Maximum number of iterations
        'MaxFunctionEvaluations', 10^10)  % Maximum number of function evaluations

    i=1;
    K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K_init, options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));

i=2;
        while i<1001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
            i=i+1;
        end

    %tau_grid=linspace(0.000001*w_w, 0.05*w_w,1000);
tau_grid=linspace(0,1500,1000);
    j=1;
    while j<1001
        tau=tau_grid(j);
        ratwtm=((w_w+tau)/w_m)^(rho/(1-rho));
        ratwtr=((w_w+tau)/r)^(theta/(1-theta));
        ratwtr2=((w_w+tau)/r)^(1/(1-theta));
        ratwtm2=((w_w+tau)/w_m)^(1/(1-rho));
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*((1/adjust*alpha_range*sigma*(l1^(1-sigma))/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            alpha_40=adjust*(40^(1-sigma))*(w_w/(sigma*l1^((1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
            alpha_range(1001:2000)=linspace(alpha_30, alpha_40, 1000);
            i=1001;
            while i<2001
            K(i)=fsolve(@(x) paramfun_lambda(x,1/adjust*alpha_range(i),sigma,theta,rho,ratrm,r,lambda,delta),K(i-1),options);
            ratK(i)=(19/K(i))^theta;
            M(i)=K(i)*(ratrm2)*l3/l4*((1-lambda)^(-1/(1-rho)))*(((1-lambda)+(lambda*ratK(i)))^e7);
            ratKW(i)=K(i)/19;
            ratKM(i)=K(i)/M(i);
            pi_rest(i)=1/adjust*alpha_range(i)*((((1-delta)*(((1-lambda)*(K(i)^theta)+lambda*(19^theta))^(rho/theta)))+(delta*(M(i)^rho)))^(sigma/rho))-(19*w_w+M(i)*w_m+r*K(i));
           i=i+1;
            end
        W_tau=((1-delta)^(sigma/(rho*(1-sigma))))*l1*((1/adjust*alpha_range*sigma/(w_w+tau)).^(1/(1-sigma)))*((l1+(l2*ratwtr))^e1)*((1+(ratwtm*l3/l4*((l1+(l2*ratwtr))^e2)))^e3);
        K_tau=(l2/l1)*ratwtr2*W_tau;
        M_tau=l3/l4*(1/l1)*ratwtm2*((l1+(l2*ratwtr))^e7)*W_tau;
        pi_tau=1/adjust*alpha_range.*((((1-delta)*(((1-lambda)*(K_tau.^theta)+lambda*(W_tau.^theta)).^(rho/theta)))+(delta*(M_tau.^rho))).^(sigma/rho))-(W_tau*(w_w+tau)+M_tau*w_m+r*K_tau);
        diff_pi=abs(pi_rest-pi_tau);
        [D2,I2]=min(diff_pi);
        
        alpha_ubar(j)=alpha_range(I2);
        if pi_rest(I2)-pi_tau(I2)<0 && I2>1
            alpha_ubar(j)=alpha_range(I2-1);
            I2=I2-1;
        end
        if I2==1000
            I2=1001;
        end
        end
        if I2==2000
            I2=1999;
        end
        
        W_next=floor((1-delta)^(sigma/(rho*(1-sigma)))*l1*(((1/adjust)*alpha_range(I2+1)*sigma/w_w)^(1/(1-sigma)))*((l1+(l2*ratwr))^((sigma-theta)/(theta*(1-sigma))))*((1+(ratwm*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((sigma-rho)/(rho*(1-sigma)))));
alpha_next=adjust*(W_next^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));
alpha_nextplus1=adjust*((W_next+1)^(1-sigma))*(w_w/(sigma*(l1^(1-sigma))))*((1-delta)^(-sigma/rho))*((l1+(l2*ratwr))^((theta-sigma)/theta))*((1+(ratwm*l3/l4*((l1+(l2*ratwr))^(rho*(theta-1)/(theta*(1-rho))))))^((rho-sigma)/rho));

        ratKW_left(j)=mean(ratKW(I2:I2)');
        ratKW_right(j)=(((1-lambda)*(w_w+tau))/(lambda*r))^(1/(1-theta));
        ratKW_discon(j)=log(ratKW_right(j))-log(ratKW_left(j));
        %Aqui es donde tengo mas dudas
        ratKM_left(j)=mean(ratKM(I2:I2)');
        ratKM_right(j)=l4/l3*l2*(ratwtr2/ratwtm2)*((l1+(l2*ratwtr))^(-e7));
        ratKM_discon(j)=log(ratKM_right(j))-log(ratKM_left(j));
        j=j+1;
    end

    %scatter(tau_grid,bunching_density);
    %scatter(tau_grid,ratKW_discon);
    %scatter(tau_grid,ratKM_discon);

    %diff_bunching=abs(bunching_density-Mccrary);
    diff_bunching=abs(ratKW_discon-RDKW);
    [D3,I3]=min(diff_bunching);
    tau_final(k,l)=tau_grid(I3);
    if D3<0.01
    ratKM_jump(k,l)=ratKM_discon(I3);
    ratKW_jump(k,l)=ratKW_discon(I3);
    bunch(k,l)=bunching_density(I3);
    else
    ratKM_jump(k,l)=9999;
    ratKW_jump(k,l)=9999;
    bunch(k,l)=9999;
    %ratKM_jump(k,l)=ratKM_discon(I3);
    %ratKW_jump(k,l)=ratKW_discon(I3);
    %bunch(k,l)=bunching_density(I3);
    end        
    clear pi_rest
    l=l+1;
    end
    k=k+1;
end

diff_KMjump=(ratKM_jump-RDKM).^2;
[D4,I4]=min(diff_KMjump);
[D5,I5]=min(D4');

rho_bygroup(iter)=rho_grid(I4(I5));
theta_bygroup(iter)=theta_grid(I5);
tau_bygroup(iter)=tau_final(I4(I5),I5);
bunching(iter)=bunch(I4(I5),I5);
KM_jump(iter)=ratKM_jump(I4(I5),I5);
KW_jump(iter)=ratKW_jump(I4(I5),I5);
end

iter=iter+1;
end
theta_median(group)=median(theta_bygroup);
theta_lb(group)=prctile(theta_bygroup,5);
theta_ub(group)=prctile(theta_bygroup,95);
rho_median(group)=median(rho_bygroup);
rho_lb(group)=prctile(rho_bygroup,5);
rho_ub(group)=prctile(rho_bygroup,95);
tau_median(group)=median(tau_bygroup);
tau_lb(group)=prctile(tau_bygroup,5);
tau_ub(group)=prctile(tau_bygroup,95);
group=group+1;
end














